from math import *
from numpy import *
import time
import random

def montecarlo_pre(funcion, limites, N):
    t0 = time.time()
    
    # Volumen
    V = prod([b - a for a,b in limites])
    # Puntos random de la funcion y funcion^2 y sus promedios
    suma_f = 0
    suma_f2 = 0
    
    for i in range(N):
        # limites = [(a, b), (c, d), ...]
        puntos = [random.uniform(a, b) for a,b in limites]
        
        f = funcion(*puntos)
        suma_f += f
        suma_f2 += f**2
        
    f_media = (1/N)*suma_f
    f2_media = (1/N)*suma_f2
    # Integral y Error
    integral = V * f_media
    error = V * sqrt((f2_media - f_media**2)/N)
    
    tf = time.time()
    
    return integral, error, (tf-t0)

def montecarlo(funcion, limites, N):
    t0 = time.time()
    
    V = prod([b - a for a,b in limites])
    dimension = len(limites)
    puntos = random.rand(N, dimension)
    
    for i,(a,b) in enumerate(limites):
        puntos[:,i] = a + (b-a)*puntos[:,i]

    f = funcion(*puntos.T)

    f_media = (1/N)*sum(f)
    f2_media = (1/N)*sum(f**2)
    
    integral = V * f_media
    error = V * sqrt((f2_media - f_media**2)/N)
    
    tf = time.time()
    
    return integral, error, (tf-t0)

# EJERCICIO 1 =================================================================
print("Ejercicio 1 " + "="*50)
print("")

def f(x):
    return  (1 - x**2)**(1.5)

lim = [(0, 1)]

print("N", "I", "Error", sep="                      ")
print("-"*55)
for n in range(2, 9):
    N = 10**n
    I, error, t = montecarlo(f, lim, N)
    print(f"{N:>12}\t{I:20.12f}\t{error:20.2e}")
print("-"*55)


# EJERCICIO 2 =================================================================
print("Ejercicio 2 " + "="*50)
print("")

lim = [(0, 1)]
def f(y):
    return (1/(y**2))*e**(-(-1 +(1/y)))

print("N", "I", "Error", sep="                      ")
print("-"*55)
for n in range(2, 9):
    N = 10**n
    I, error, t = montecarlo(f, lim, N)
    print(f"{N:>12}\t{I:20.12f}\t{error:20.2e}{t:10.2f}s")
print("-"*55)


# EJERCICIO 3 =================================================================
print("Ejercicio 3 " + "="*50)
print("")

def funcion_radio(x, y, R):
    if x**2 + y**2 <= R**2:
        return 1
    else:
        return 0

def area_circulo(R, N):
    t0 = time.time()
    
    x = random.uniform(-R,R,N)
    y = random.uniform(-R,R,N)
    suma = 0
    suma2 = 0
    for i in range(N):
        r = funcion_radio(x[i], y[i], R)
        suma += r
        suma2 += r**2
    
    tf = time.time()
    area = (4*R**2)/N *suma
    error = (4*R**2)*sqrt(((suma2/N)-(suma/N)**2)/N)
    return area, error, (tf-t0)

R = 1.5
print("N", "I", "Error", sep="                      ")
print("-"*55)
for n in range(2, 9):
    N = 10**n
    I, error, t = area_circulo(R, N)
    print(f"{N:>12}\t{I:20.12f}\t{error:20.2e}{t:10.2f}s")
print("-"*55)


# EJERCICIO 4 =================================================================
print("Ejercicio 4 " + "="*50)
print("")

R = 1.66

def funcion_esfera(x, y, z, R):
    if x**2 + y**2 + z**2 <= R**2:
        return 1
    else:
        return 0

def volumen_esfera(R, N):
    
    t0 = time.time()

    x = random.uniform(-R, R, N)
    y = random.uniform(-R, R, N)
    z = random.uniform(-R, R, N)
    suma = 0
    suma2 = 0

    for i in range(N):
        r = funcion_esfera(x[i], y[i], z[i], R)
        suma += r
        suma2 += r**2

    Vcubo = (2*R)**3
    volumen = Vcubo * suma/N
    error = Vcubo * sqrt((suma2/N - (suma/N)**2)/N)

    tf = time.time()

    return volumen, error, (tf-t0)

print("N", "I", "Error", sep="                      ")
print("-"*55)
for n in range(2, 9):
    N = 10**n
    V, error, t = volumen_esfera(R, N)
    print(f"{N:>12}\t{V:20.12f}\t{error:20.2e}{t:10.2f}s")
print("-"*55)


# EJERCICIO 5 =================================================================
print("Ejercicio 5 " + "="*50)
print("")

R, r = 3, 1
limites = [(1, 4), (-3, 4), (-1, 1)]

def toro(x, y, z):
    if (sqrt(x**2 + y**2) - R)**2 + z**2 <= r**2:
        return 1
    else:
        return 0



print("N", "I", "Error", sep="                      ")
print("-"*55)
for n in range(2, 9):
    N = 10**n
    I, error, t = montecarlo_pre(toro, limites, N)
    print(f"{N:>12}\t{I:20.12f}\t{error:20.2e}\t{t:10.2f}s")
print("-"*55)